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Abstract 

We propose a model of ripples in suspended graphene sheets based on plate equations that are 
made discrete with the periodicity of the honeycomb lattice and then periodized. In addition, 
the equation for the displacements with respect to the planar configuration contains a double-well 
site potential, a nonlinear friction and a multiplicative white noise term satisfying the fluctuation- 
dissipation theorem. The nonlinear friction terms agree with those proposed by Eichler et al 
[Nature Nanotech. 6, 339 (2011)] to explain their experiments with a graphene resonator. The site 
double-well potential indicates that the carbon atoms at each lattice point have equal probability 
to move upward or downward off-plane. For the considered parameter values, the relaxation time 
due to friction is much larger than the periods of membrane vibrations and the noise is quite 
small. Then ripples with no preferred orientation appear as long-lived metastable states for any 
temperature. Numerical solutions confirm this picture. 

PACS numbers: 61.48. Gh,68.65.Pq,05.45.-a 
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I. INTRODUCTION 



The first visualizations of atoms in suspended graphene sheets showed that they were 
covered with ripples.- These ripples are several nanometers long waves of the sheet without 
a preferred direct ion. The initial preparation of a suspended graphene bilayer may result in 
Angstrom-sized ripples that are absent in graphene monolayers.- In experiments on graphene 
sheets suspended on substrate trenches, there appear much longer and taller waves (close 
to a micron scale) directed parallel to the applied stress. 4 These long wrinkles are thermally 
induced and can be explained by continuum elasticity- Ripples are expected to be important 
for electronic transport in graphene,- and there is active research about the effects of ripples 
and strain on electronic properties, including possible strain engineering 

The earliest theoretical studies of ripples using Monte Carlo^ or molecular dynamics 
simulations^ have shown that ripples may be connected to variable length a bonds of carbon 
atoms and may be due to thermal fluctuations. Other studies explored the connection 
between rippling and electronic properties" and suggested that, at zero temperature, 
the electron-phonon coupling may drive the graphene sheet into a quantum critical point 
characterized by the vanishing of the bending rigidity of the membrane.— There is a buckling 
transition between the regimes of positive and negative tension in graphene. Then the 
inherent tension in suspended graphene produces ripples with a characteristic length and 
aspect ratio.— Thus rippling already occurs at OK and it is due to the strong interplay 
between the dynamics of the free electrons in the membrane and its mesoscopic structure. 
While most of rippling studies are focused in equilibrium behavior, in this paper we want to 
analyze the dynamics of ripples using a model based on discrete elasticity that incorporates 
nonlinear friction and white noise sources. 

When observing ripples with a high-resolution transmission electron microscope, the 
graphene sample is bombarded by a low-intensity electron beam. Even though the elec- 
tron beam does not have sufficient energy to knock off carbon atoms, it may create defects 
by inducing bond rotation and certainly excite the atoms. Thus an observed graphene 
sample is continuously excited and cannot be considered to be in equilibrium.-^ We have 
suggested that the electron beam may push the carbon atoms vertically away from the pla- 
nar configuration of the graphene sheet in a random fashion and its effect could be modeled 
by coupling each carbon atom to an Ising spin that tries to move it off-plane.— In a simpler 
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context, one- dimensional (ID) mechanical systems coupled to Ising spins that flip randomly 
according to Glauber dynamics exhibit ripples at any temperature in the limit as the spin 
relaxation time is much longer than the vibration periods of the mechanical system.— In 
this spin-string system, the ripples are metastable quasi-equilibrium states about which the 
mechanical system experiences rapid oscillations. Both the ID spin-string and the spin- 
membrane (graphene) models exhibit a buckling transition at a critical temperature, above 
which the string or the membrane is flat. Below the critical temperature, the spins are 
polarized and the mechanical part of the system (string or membrane) is buckled. 

In the case of suspended graphene sheets, there is a large degree of arbitrariness in the 
assignation of numerical values to the spin-atom coupling constant and to the spin-flip rate^ 
and there is not yet unambiguous experimental evidence of a buckling transition at finite 
temperature. Instead, quantum theories point out that rippling due to the strong inter- 
play between the free electron dynamics in the membrane and its mesoscopic structure,— 1 ^ 
already occurs at OK.— Thus it is convenient to have a model exhibiting rippling at any 
temperature and whose parameters can be calibrated with experiments. In this paper we 
present such a model. 

We describe the position of carbon atoms by discrete elasticity equations previously used 
to model defect dynamics in graphene.— - — Neighboring carbon atoms in a graphene sheet 
attach to each other using three of their bonds. The fourth bond is not saturated and, similar 
to the effect of the electron beam, it may try to pull the atom upward or downward from the 
flat sheet configuration. This trend may be modeled by placing each atom in a double-well 
potential. This is reminiscent of the system of a soft spin coupled to elasticity used to study 
two-dimensional (2D) melting mediated by dislocations. 19 Without the double- well potential, 
our numerical simulations of the time-dependent von Karman plate equations (discretized 
on the honeycomb lattice) show that a suspended sheet that is initially in a rippled state 
remains rippled whereas a flat sheet remains flat. In the absence of longitudinal stretching, 
the initial conditions of the plate equations select whether the graphene sheet remains flat 
or it acquires curvature or ripples. In experiments, rippling is always observed, so we have 
to look for a mechanism to produce it spontaneously. Theories including electron-phonon 
coupling indicate that, for infinitesimal stress, ripples are spontaneously formed even at 
OK,— which introduces an extra length similar to that in our double-well potential. In 
addition, we add nonlinear friction and an appropriate white noise forcing that satisfies the 
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fluctuation-dissipation theorem. The fundamental mode of the equations of motion can be 
projected on a Duffing oscillator with nonlinear friction and white noise forcing. It tuns 
out that Eichler et al have interpreted their measurements of graphene resonators by using 
such a Duffing model (without white noise forcing). 20 In this latter work, measurements 
of carbon nanotube and graphene resonators are used to calculate the coefficients of the 
Duffing oscillator. Assuming that the minima of the potential are upward or downward 
a distance wq off-plane, we use their coefficients to fit the parameters in the equations of 
motion of our graphene model. The only free parameter in our model is Wq. While a more 
elaborated theory could derive a temperature and strain dependent wq from first principles 
(for example, considering an electron-phonon interaction)™, we have considered fixed values 
of Wo from 0.5 to 1 nm for the sake of simplicity. Quantum theories produce lengths of this 
order.— ^ Numerical solutions of our model yield ripples without preferred orientation in 
agreement with experimental observations.- 1 ^ 

The rest of the paper is as follows. The model we use is described in Section [Til In Section 
IHIl we calculate numerical values of the coefficients used in our model by fitting it to one of 
the graphene resonators studied by Eichler et al.— In Section IIVI we discretize the equations 
of motion on the honeycomb lattice, write them in primitive coordinates and periodize them, 
so that integer multiples of displacements on primitive directions leave the lattice unchanged. 
We then solve numerically the resulting stochastic equations of motion. We show that long- 
lived ripples with no preferred orientation emerge after a short transient even from random 
initial conditions. These ripples correspond to those observed in experiments.- 1 ^ Lastly, 
Section [V] contains our conclusions. 

II. MODEL 

In the graphene sheet, carbon atoms have o bond orbitals constructed from sp 2 hybrid 
states oriented in the direction of the bond that accommodate three electrons per atom. 
The other electrons go to p states oriented perpendicularly to the sheet. These orbitals bind 
covalently with neighboring atoms and form a narrow ir band that is half-filled. The presence 
of bending and ripples in graphene modifies its electronic structure.— Out-of-plane convex 
or concave deformations of the sheet have in principle equal probability and transitions 
between these deformations are associated with the bending energy of the sheet. A simple 
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way to model this situation is to consider that out-of-plane deformations are described by 
a double-well site potential that tries to set vertical deflections of the sheet, w(x,y) to ±u> 
and contributes the free energy: 
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where p 2 is the 2D mass density (mass per unit area) and (p has units of velocity square. The 
elastic free energy of the graphene sheet in the continuum limit is that of a 2D membrane,— 

F a = \ I fi(V 2 w) 2 + {Xul + 2~pu 2 k )} dxdy, 



(2) 

Uik = ^(d Xk Ui + d Xi u k + d Xi wd Xk w), i,k = l, 2, (3) 

where (ui,u 2 ) = (u(x,y),v(x,y)), and k, A and p are the in-plane displacement vector, the 
bending stiffness (measured in units of energy) and the 2D Lame coefficients of graphene 
(measured in units of force per unit length), respectively. V = (d x , d y ) is the 2D gradient and 
V 2 the 2D laplacian. In we have ignored the small in-plane nonlinear terms d Xi ud Xk u + 
d Xi vd Xk v. 

From the total free energy F = F„ + F DW , we obtain the equations of motion 
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(£(x, y, t)) = 0, (£(x, y, t)£{x, y, t)) = 5(x - x')5{y - y')5(t - t'), (7) 

where P is the membrane stress, 9 is the temperature measured in units of energy and — (7 + 
fjw 2 )d t w is a nonlinear friction force used by Eichler et al 20 to interpret their experiments 



with a forced damped graphene resonator. The intensity y26 l (7 + fjw 2 ) of the white noise 
is related to the friction by the fluctuation-dissipation theorem. All the parameters A, 
p, p 2 , k, P, u)q, (p, 7, 77 and 6 are positive. 



III. PARAMETER IDENTIFICATION 



To identify the parameters A, jl, p 2 , k, P, wq, ip x 7, rj and 9, we shall consider the 



graphene resonator corresponding to Figure 3 in Ref. 
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The temperature is 4K, so that we 
use A = 3.25 eV/A 2 and jl = 9.44 eV/A 2 corresponding to OK in Table 1 of Ref. I23J. The 
graphene resonator has length L = 1.7 pm, width W = 120 nm and is pinned at x = and 
x = L. The sides at y = and y = W are free.— Let us assume that u — v — 0, ignore 
bending and noise, and assume that w is in the lowest possible eigenstate: 

7TX 

w(x, y, t) = sin — . (8) 

We now insert Eq. (JSJ) in (EJ) , multiply the result by sin(7ra;/L), integrate and divide by 
J Q L sin 2 (773; /L)dx = L/2. We obtain 

The stress P in the sheet is due to the inherent stretching of a stressed planar membrane, 
P = E2S0 [where Eq is the homogeneous strain and E 2 = 4/2(A+/i) / (A+2/i) is the 2D Young's 
modulus], and an additional stress AP due to its bending. The latter is proportional to the 
variation of the membrane area: 

AS f L f W Jl+ \WW\ 2 - 1 , , 1 f L f W ^ , 2J , tt 2 

in which Eq. (JS} has been used. Thus the total stress P is 
which inserted in Eq. ([9]) produces 



rft 2 \p 2 L 2 w 2 \Awl 8p 2 L 4K ^' 4p 2 L 4 

+ (12) 



p 2 Ap 2 J dt 

This should be identified with the nonlinearly damped Duffing oscillator used by Eichler 
et^Si as a model of their graphene resonator: 

k a 3 7* + V*x 2 . , . 

x = x x x, (13) 

mm m 



FIG. 1: (Color online) Neighbors of a given atom A in sublattice 1 (dark blue). The primitive 
vectors are a, b and c that connect atom B to three different next nearest neighbors in its same 
sublattice. 



where k/m = 7r 2 To/(mL) is the resonant frequency squared, m = 3.9 x 10 _19 kg is the mass 
of the sheet, T = 110 nN is the intrinsic tension, a = 1.4 x 10 16 kg m~ 2 s~ 2 , 7* = 8.7 x 1CT 14 
kg/s, and rf = 1.5 x 10 7 kg m~ 2 s -1 . We find 



9 = 4a _ AX + 2g 1 ^(A + A) = 4 3142 x w u m _ 2s _ 2 (M) 
ti>3 3m p 2 L' \2 3(A + 2/i) 2 ' 



Setting wq = 5 A, we obtain = 2.6964 x 10 3 m 2 s 2 , an intrinsic strain Eo = 0.0010067 
and stress Pq = 0.3492 kg s~ 2 . The friction coefficients are 7//O2 = 2.2308 x 10 5 s _1 and 
f}jp 2 = 5.1282 x 10 25 m- 2 s- 1 . 



IV. NUMERICAL RESULTS 
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we discretize the equations of 



Following the same procedure as in References 
motion on the honeycomb lattice of graphene. Atoms belonging to the two sublattices 
comprising the lattice are distinguished in Figure [TJvith light and dark colors, although they 



are all identical carbon atoms. The three nearest and six next nearest neighbors of the atom 
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A with coordinates (x, y) in Fig. [TJ are 



a a \ fa a 

n i = \ x -niV - h n 2= \x + -,y- 



2y/3f Z V 2 2V3/ 

o \ ( a a\/3 
n 3 = ( x,y + -j= I, ra 4 =lx--,y — 

n 5 =|x + -,y —I, 7i 6 = (ar-o,j/), 

, a a\/3 
n 7 =(x + a,y), n 8 = \ x - ~,y + 



2 7J 2 



, a a\/3\ 

% = I x + 2'^ + )• ( 15 ) 

To discretize the equations of motion using the honeycomb lattice, we define the following 
operators acting on functions of the node A coordinates (x,y), 

Tu = [u(n x ) - u{A)\ + [u(n 2 ) - u{A)\ 

+ [u(n 3 )-u(A)]~(dlu + d*u)^, (16) 

Hu = [u{n 6 ) - u{A)\ + [u{n 7 ) - u{A)\ ~ a 2 d 2 x u, (17) 

D\U = \u{n A ) - u{A)\ + [u(n 9 ) - u(A)] (18) 

1 -.2.. , V3 o o „, , 3 o 1 2 



- c^w + — d x d y u + - c^it j (/ 
D 2 u = [u(n 5 ) - u(A)} + [u(n 8 ) - u(A)} (19) 

A/jti = 14(717) — U (A) ~ (c^tt) a, (20) 

A,,m = u(n 3 ) - u(A) ~ (d y u) -^=, (21) 

Bw = [Twfa) - Tw(A)] + [Tw{n 2 ) - Tw(A)} 

+ [Tw(n 3 ) - T W (A)] ~ ^(«9 2 + d 2 ) 2 w. (22) 

with similar definitions for points B in sublattice 2. Taylor expansions of these finite dif- 
ference combinations about (x, y) yield the partial derivative expressions written above as 
a —7- 0. We now replace the derivatives in (j4]), (jSJ) and (jSJ), <9^w, d 2 u, d x d y u, d x w, d y w and 
(V 2 ) 2 w by the operators #w/a 2 , (AT - H)u/a 2 , (D 1 - D 2 )u/(V3a 2 ), A h w/a, y/3A v w/a 
and 16-Biy/a 4 , respectively, with similar substitutions for the derivatives of v and w. The 
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resulting equations at each point of the lattice are 



p 2 a 2 d 2 t u = AjlTu + (A + fi) Hu + ^£ (D 1 - D 2 )v 
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(25) 



Here 



(26) 



is the stress in the membrane (see section IHIj) and iV = J2 X y 1 * s the total number of atoms 
in the graphene sheet. In the limit a — > 0, the continuous delta function 5(x — x') becomes 
5 xx '/a, where 8 XX > = 1 if x = x' and otherwise. Then the discretized white noise in Eq. 
(I25p . a£(x,y,t), has zero mean and correlation 



(a£(x, y, t) a£(x', y', t')) = 5 xx >5 yy i5(t - t') 



(27) 
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Possible defects inserted in the graphene sheet are the cores of dislocations. To account for 
them, we have to write these equations of motion in primitive coordinates and periodize all 
difference operators appearing in them along primitive directions. The resulting equations 
of motion become those in Ref. Q for w = 0. If there are not defects in the graphene sheet, 
the equations of motion (T^3|)-([^3J) or their periodized version produce the same results. 



u, v, w, Wq t £ 8 P 



\+2fi 



,,,,,1/4 

^» (A + 2p) a 2 A + 2/i 



A 10- 14 s m-h- 1 / 2 eV eV/A 2 

1.42V3 1.11109 3.85723 x 10 16 133.8688 22.13 

TABLE I: Units used to nondimensionalize equations. 

It is convenient to write the previous equations of motion in nondimensional units. Using 
Table HJ we obtain the following nondimentional equations of motion 

d 2 u = A^Tu + (A + fi) (hu + (Dl ~J? 2)V ^J 

+ (A + /i) [A h w Hw + A v w (Di — D 2 )w} + 4/x A h w Tw, (28) 
d 2 v = ATv + (A + fi) ( ^ Dl 'J^^ - Hv j + 4^ A v w Tw 

+ h^£[A h w{D 1 - D 2 )w - SAyW Hw], (29) 
v3 

d 2 t w = [Hu + 2A v w (D 1 - D 2 )w + A h w Hw] A h w 



+ 



V3(4T - H)v + 3A v w (AT - H) 



w 



A v w 



+ (A + fi)(D 1 - D 2 )u A v w + —= A h w (D 1 - D 2 )v 

a/3 

+ n A v w [V3(AT - H)u + VSHv + (D 1 - D 2 )v] + 2A h u (2AT + fxH) 

+ 2v / 3A^ [2 AT + /i(4T - H)]w + 2fi (^A v u + ^ J (D 1 - D 2 )w 

+ [( A h w) 2 + (A v w) 2 ] (2AT + (iH)w + A/j Tw (AT - H)w + AP Tw 

1 o)—? — h + r]w 2 )d t w+ \/29h + -qw 2 ) E(x, y, t). 

w 2 J w 2 



(30) 
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where u, v, w and t are now nondimensional and 



p = frS ( e ° + 2* ^ l(A " 1 " )2 + 9(A " W)2] I ' (31) 



X » = rAz, (32) 



A + 2/i' ' A + 2/i 



(S(x,y,*)) = 0, (S(x,y,t)S(x',y',t')> = <W<W*(f-0- (33) 
The dimensionless parameters appearing in these equations are listed in Table [Til We have 



used the values of the 2D Lame moduli at OK given in Ref. |23|: A + 2/2 = 22.13 eV/A 2 , 
p = 9.44 eV/A 2 . The bending rigidity k = 0.8 eV has been extrapolated from Fig. 7 of 
Ref. [24j. The 2D mass density is p2 = 7.236 x 10 -7 kg/m 2 , and the parameters identified in 
Section [TTTJ are <p/w^ = 4.3142 x 10 34 hi'V 2 and, for w = 5 A, we have e = 0.0010067, 
P Q = 0.3492 kg s" 2 , 7/p 2 = 2.2308 x 10 5 s^ 1 and fj/p 2 = 5.1282 x 10 25 m- 2 s" 1 . In our 
simulations, we have considered shorter samples than in Eichler et al's experiments for the 
sake of computational expediency. To keep the same frequency, we need to rescale (p as L~ 4 , 
and both rj and 7 as 1/L. For a 80 x 80 lattice, this yields the parameters in the second row 
of Table HQ 



wo 8 


k 7 r] 




£0 


vio 6 


16k 7a V a3 






a (A+2/i)a 2 


(X+2p,)a 2 ^/ P2 (A+2/i) VP2(A+2A) 






2.03292 2.57484 x 10~ 6 


0.095616 0.00342536 0.0476338 


5.5028 x 10" 12 


0.0010067 


2.03292 2.57484 x 10" 6 


0.095616 0.295949 4.11553 


3.06636 x 10~ 4 


0.0010067 



TABLE II: Dimensionless parameters for Fig. 3 of Ref. [2CJ with wq = 5 A. The numbers in the last 



row correspond to a shorter computational lattice with L = 80a. 

The model equations produce spontaneous ripple generation in which ripples have no 
preferred orientation and their size is comparable with that observed in experiments. We 
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have considered that the graphene sheet has a strain similar to that in Ref. 
subjacent strain flattens the ripples as compared to those in an unstrained sheet. Initially 
the graphene sample is flat. As shown in Fig. E^a), random forcing pulls atoms off-plane 
and randomly oriented ripples that are very small in size and height form. After a short 
transient stage, the ripples have reached a quasi-stationary state in which they are domains 
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FIG. 2: (Color online) Ripples formed in a suspended 80 x 80-honeycomb graphene sheet (with 
12,800 atoms), (a) Density plot of the off-plane displacement as a function of in-plane coordinates 
at t = 1. (b) The same at t = 50, and (c) t = 100. (d) Average density plot over all times from 
t = 1 to 100 with unit time step. Parameter values as in the second row of Table ITTl 

of atoms displaced variable distances either upward or downward the horizontal plane, cf. 
Fig. E](b) and (c) and Fig. [3] (the latter, for a 160 x 160 lattice). Ripple domains change 
slowly in height and size, essentially due to displacements of the atoms in their boundaries. 
This is seen when plotting the time averaged density plot of Fig. [2]^d). More clearly, the 
averages between more advanced times (between 300 and 400) do not change much from the 
instantaneous time resolved density plots, as shown in Fig. H]for a larger 300 x 300 graphene 
lattice. 

In the quasi-stationary state, ripple size and height increase with sample size and with 
Wq. For instance, in a 80 x 80 lattice, ripples are 2.5 nm long and 0.1 pm high for wq = 5A. 
Increasing w a factor 10 does not change the ripple height that much but there appear long 
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FIG. 3: (Color online) Ripples in the central part of a 160 x 160 graphene lattice (with 51,200 
atoms) . 

connected domains with atoms either up or down the horizontal plane. For longer 160 x 160 
samples with wq = lnm as in Fig. [31 the ripple length has increased to about 10 nm with 
height of about 0.5 pm. Increasing sample size to 300 x 300 results in longer ripples (15 nm) 
and height about 1 pm, as shown in Fig. |H As graphene sheets are suspended on holes that 
have 1 fim diameter, the effects of boundaries are felt comparatively less by atoms far from 
them and the ripple heights increase correspondingly to several Angstroms.-^ 

V. CONCLUSIONS 

We have proposed a mechanism to spontaneously produce ripples in a suspended graphene 
sheet. Besides the membrane free energy, carbon atoms are placed in a double- well potential 
at each lattice site. Thus they experience a force that tries to displace them vertically off- 
plane, either up or down with equal probability. Inspired by experiments with graphene 
resonators, we assume that the atoms undergo a nonlinear friction force and are subject 
to the corresponding multiplicative white noise force satisfying the fluctuation-dissipation 
theorem. The resulting time scale of friction is much larger than the vibration periods of 
the membrane. Both the site double-well potential and the site nonlinear friction should 
be considered as simple mechanisms that provide a nonlinearly damped Duffing oscillator 
equation for the amplitude of the leading membrane normal mode. Extra soft-spin fields^ 
or a viscous component of the stress tensor plus a white noise source as in fluctuating 
hydrodynamics^ may prove better options and yet yield a similar Duffing oscillator equation. 
We have solved numerically the equations of the model using parameters obtained from ab- 
initio calculations and from experiments in the literature on graphene resonators. Even 
from initial conditions corresponding to a flat membrane, the nonlinear force created by 
the double-well potential and the white noise forcing produce stable rippling after a short 
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FIG. 4: (Color online) Same as in Fig. Q] for a 200 x 200 honeycomb lattice with wq = 1 nm 
(180,000 atoms). Times are (a) t = 10, (b) 250, (c) 510. (d) Average density plot from t = 410 to 
510, using a t = 10 time step. Parameters are as in Table U except for the rescaled 7 = 0.11838, 
77 = 1.64621, 99 = 1.25598 x 10~ 4 . 

transient. The ripples are formed by domains all whose atoms are displaced upward or 
downward off plane. Once they have acquired a sufficient size, of the order of experimentally 
observed ripples, the domains vary slowly by annexing or losing atoms in their periphery. In 
this scenario, randomly oriented nanometer-sized ripples with no preferred direction appear 
as long-lived metastable states for any temperature. 
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